//-----------------------------------------------------------------------------
// Sun.java
//
// (c) 2004 Jonathan Stott
//
// Created on 30-Mar-2004
//
// 0.3 - 06 Mar 2005
//  - Updated incorrect sign of longitude (should be negative for west)
// 0.2 - 13 Apr 2004
//  - Change time zones to use the TimeZone class
//  - Adjust hours greater than 23 in convertTime()
// 0.1 - 30 Mar 2004
//  - First version 
//-----------------------------------------------------------------------------

package uk.me.jstott.sun;

import java.util.TimeZone;

import uk.me.jstott.coordconv.LatitudeLongitude;


/**
 * Provides methods to calculate sunrise, sunset and civil, nautical and
 * astronomical twilight times for a given latitude and longitude.
 *
 * Based on functions available at
 * http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html
 *
 * For more information on using this class, look at
 * http://www.jstott.me.uk/jsuntimes/
 * 
 * @author Jonathan Stott
 * @version 0.2
 */
public class Sun {
  
  private static final double SUNRISE_SUNSET_ZENITH_DISTANCE = 90.83333;
  private static final double CIVIL_TWILIGHT_ZENITH_DISTANCE = 96.0;
  private static final double NAUTICAL_TWILIGHT_ZENTITH_DISTANCE = 102.0;
  private static final double ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE = 108.0;


  /**
   * Calculate the UTC of a morning phenomenon for the given day at the given
   * latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param latitude latitude of observer in degrees
   * @param longitude longitude of observer in degrees
   * @param zenithDistance one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE,
   *                       Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE,
   *                       Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE,
   *                       Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE.
   * @return time in minutes from zero Z
   */
  private static double morningPhenomenon(double julian,
                                          double latitude,
                                          double longitude,
                                          double zenithDistance) {
   // longitude = -longitude;
    double t = julianDayToJulianCenturies(julian);
    double eqtime = equationOfTime(t);
    double solarDec = sunDeclination(t);
    double hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
    double delta = longitude - deg(hourangle);
    double timeDiff = 4 * delta;
    double timeUTC = 720 + timeDiff - eqtime;

    // Second pass includes fractional julian day in gamma calc
    double newt = julianDayToJulianCenturies(
        julianCenturiesToJulianDay(t) + timeUTC / 1440);
    eqtime = equationOfTime(newt);
    solarDec = sunDeclination(newt);
    hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
    delta = longitude - deg(hourangle);
    timeDiff = 4 * delta;
    
    return 720 + timeDiff - eqtime;
  }
  
  
  /**
   * Calculate the UTC of sunrise for the given day at the given latitude and
   * longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use (e.g. Sun.GMT or Sun.PST)
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time sunriseTime(double julian,
                                   LatitudeLongitude ll,
                                   TimeZone timeZone,
                                   boolean dst) {
    double timeMins = morningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        SUNRISE_SUNSET_ZENITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Calculate the UTC of morning civil twilight for the given day at the given
   * latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use (e.g. Sun.GMT or Sun.PST)
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time morningCivilTwilightTime(double julian,
                                              LatitudeLongitude ll,
                                              TimeZone timeZone,
                                              boolean dst) {
    double timeMins = morningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        CIVIL_TWILIGHT_ZENITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Calculate the UTC of morning nautical twilight for the given day at the
   * nautical latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use (e.g. Sun.GMT or Sun.PST)
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time morningNauticalTwilightTime(double julian,
                                                 LatitudeLongitude ll,
                                                 TimeZone timeZone,
                                                 boolean dst) {
    double timeMins = morningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        NAUTICAL_TWILIGHT_ZENTITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Calculate the UTC of morning astronomical twilight for the given day at the
   * given latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use (e.g. Sun.GMT or Sun.PST)
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time morningAstronomicalTwilightTime(double julian,
                                                     LatitudeLongitude ll,
                                                     TimeZone timeZone,
                                                     boolean dst) {
    double timeMins = morningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Calculate the UTC of an evening phenomenon for the given day at the given
   * latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param latitude latitude of observer in degrees
   * @param longitude longitude of observer in degrees
   * @param zenithDistance one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE,
   *                       Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE,
   *                       Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE,
   *                       Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE.
   * @return time in minutes from zero Z
   */
  private static double eveningPhenomenon(double julian,
                                          double latitude,
                                          double longitude,
                                          double zenithDistance) {        
    double t = julianDayToJulianCenturies(julian);

    // First calculates sunrise and approx length of day
    double eqtime = equationOfTime(t);
    double solarDec = sunDeclination(t);
    double hourangle = hourAngleEvening(latitude, solarDec, zenithDistance);

    double delta = longitude - deg(hourangle);
    double timeDiff = 4 * delta;
    double timeUTC = 720 + timeDiff - eqtime;

    // first pass used to include fractional day in gamma calc
    double newt = julianDayToJulianCenturies(
        julianCenturiesToJulianDay(t) + timeUTC / 1440);
    eqtime = equationOfTime(newt);
    solarDec = sunDeclination(newt);
    hourangle = hourAngleEvening(latitude, solarDec, zenithDistance);

    delta = longitude - deg(hourangle);
    timeDiff = 4 * delta;
    
    return 720 + timeDiff - eqtime;
  }
  
  
  /**
   * Calculate the UTC of sunset for the given day at the given latitude and
   * longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use (e.g. Sun.GMT or Sun.PST)
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time sunsetTime(double julian,
                                  LatitudeLongitude ll,
                                  TimeZone timeZone,
                                  boolean dst) {
    double timeMins = eveningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        SUNRISE_SUNSET_ZENITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Calculate the UTC of evening civil twilight for the given day at the given
   * latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time eveningCivilTwilightTime(double julian,
                                              LatitudeLongitude ll,
                                              TimeZone timeZone,
                                              boolean dst) {
    double timeMins = eveningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        CIVIL_TWILIGHT_ZENITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Calculate the UTC of evening nautical twilight for the given day at the
   * given latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use (e.g. Sun.GMT or Sun.PST)
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time eveningNauticalTwilightTime(double julian,
                                                 LatitudeLongitude ll,
                                                 TimeZone timeZone,
                                                 boolean dst) {
    double timeMins = eveningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        NAUTICAL_TWILIGHT_ZENTITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Calculate the UTC of evening astronomical twilight for the given day at the
   * given latitude and longitude on Earth
   *
   * @param julian Julian day
   * @param ll latitude and longitude of observer in degrees
   * @param timeZone the timeZone to use (e.g. Sun.GMT or Sun.PST)
   * @param dst true if daylight savings time (summer time) should be taken into
   *            account
   * @return time in minutes from zero Z
   */
  public static Time eveningAstronomicalTwilightTime(double julian,
                                                     LatitudeLongitude ll,
                                                     TimeZone timeZone,
                                                     boolean dst) {
    double timeMins = eveningPhenomenon(julian,
                                        ll.getLatitude(),
                                        -ll.getLongitude(),
                                        ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE) +
                      (timeZone.getRawOffset() / 60000.0);
    if (dst) timeMins += 60.0;
    
    return convertTime(timeMins);
  }
  
  
  /**
   * Convert a double representing a time to a Time object
   *
   * @param time
   * @return
   */
  private static Time convertTime(double time) {
    int hours = (int)(time / 60.0);
    int minutes = (int)(time - (hours * 60));
    double seconds = (time - minutes - (hours * 60)) * 60;
    if (hours > 23) hours %= 24;
    
    //Nekloth's correction
    int seconds_int;
    if (hours<0) hours += 24;
    if (minutes<0) minutes = (int) minutes + 60;
    if (seconds<0) {
    	seconds_int = (int) (seconds+60);
    } else {
    	seconds_int = (int) seconds;
    }
    return new Time((int)hours, (int)minutes, seconds_int);
  }
  
  
  /**
   * Convert Julian Day to centuries since J2000.0
   * 
   * @param julian The Julian Day to convert
   * @return the value corresponding to the Julian Day
   */
  private static double julianDayToJulianCenturies(double julian) {
    return (julian - 2451545) / 36525;
  }
  
  
  /**
   * Convert centuries since J2000.0 to Julian Day
   * 
   * @param t Number of Julian centuries since J2000.0
   * @return The Julian Day corresponding to the value of t
   */
  private static double julianCenturiesToJulianDay(double t) {
    return (t * 36525) + 2451545;
  }
  
  
  /**
   * Calculate the difference between true solar time and mean solar time
   * 
   * @param t Number of Julian centuries since J2000.0
   * @return 
   */
  private static double equationOfTime(double t) {
    double epsilon = obliquityCorrection(t);
    double l0 = geomMeanLongSun(t);
    double e = eccentricityOfEarthsOrbit(t);
    double m = geometricMeanAnomalyOfSun(t);
    double y = Math.pow((tan(rad(epsilon) / 2)), 2);
    
    double Etime = y * sin(2 * rad(l0)) - 2 * e * sin(rad(m)) + 4 * e * y *
                   sin(rad(m)) * cos(2 * rad(l0)) - 0.5 * y *
                   y * sin(4 * rad(l0)) - 1.25 * e * e * sin(2 * rad(m));

    return Math.toDegrees(Etime) * 4;      
  }
  
  
  /**
   * Calculate the declination of the sun
   * 
   * @param t Number of Julian centuries since J2000.0
   * @return The Sun's declination in degrees
   */
  private static double sunDeclination(double t) {
    double e = obliquityCorrection(t);
    double lambda = sunsApparentLongitude(t);

    double sint = sin(rad(e)) * sin(rad(lambda));
    return deg(asin(sint));
  }
  
  
  /**
   * calculate the hour angle of the sun for a morning phenomenon for the given
   * latitude
   * 
   * @param lat Latitude of the observer in degrees
   * @param solarDec declination of the sun in degrees
   * @param zenithDistance zenith distance of the sun in degrees
   * @return hour angle of sunrise in radians
   */
  private static double hourAngleMorning(double lat,
                                         double solarDec,
                                         double zenithDistance) {
    return (acos(cos(rad(zenithDistance))
        / (cos(rad(lat)) * cos(rad(solarDec))) - tan(rad(lat))
        * tan(rad(solarDec))));
  }
  
  
  /**
   * Calculate the hour angle of the sun for an evening phenomenon for the given
   * latitude
   *
   * @param lat Latitude of the observer in degrees
   * @param solarDec declination of the Sun in degrees
   * @param zenithDistance zenith distance of the sun in degrees
   * @return hour angle of sunset in radians
   */
  private static double hourAngleEvening(double lat,
                                         double solarDec,
                                         double zenithDistance) {
    return -hourAngleMorning(lat, solarDec, zenithDistance);        
  }
  
  
  /**
   * Calculate the corrected obliquity of the ecliptic
   * 
   * @param t Number of Julian centuries since J2000.0
   * @return Corrected obliquity in degrees
   */
  private static double obliquityCorrection(double t) {
    return meanObliquityOfEcliptic(t) + 0.00256
        * cos(rad(125.04 - 1934.136 * t));          
  }
  
  
  /**
   * Calculate the mean obliquity of the ecliptic
   *
   * @param t Number of Julian centuries since J2000.0
   * @return Mean obliquity in degrees
   */
  private static double meanObliquityOfEcliptic(double t) {
    return 23 + (26 + (21.448 - t * (46.815 + t * (0.00059 - t * (0.001813)))
        / 60)) / 60;
  }
  
  
  /**
   * Calculate the geometric mean longitude of the sun
   * 
   * @param number of Julian centuries since J2000.0
   * @return the geometric mean longitude of the sun in degrees
   */
  private static double geomMeanLongSun(double t) {
    double l0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
    
    while ((l0 >= 0) && (l0 <= 360)) {
      if (l0 > 360) {
        l0 = l0 - 360;
      }
      
      if (l0 < 0) {
        l0 = l0 + 360;
      }
    }
          
    return l0;
  }
  
  
  /**
   * Calculate the eccentricity of Earth's orbit
   * 
   * @param t Number of Julian centuries since J2000.0
   * @return the eccentricity
   */
  private static double eccentricityOfEarthsOrbit(double t) {
    return 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
  }
  
  
  /**
   * Calculate the geometric mean anomaly of the Sun
   *
   * @param t Number of Julian centuries since J2000.0
   * @return the geometric mean anomaly of the Sun in degrees
   */
  private static double geometricMeanAnomalyOfSun(double t) {
     return 357.52911 + t * (35999.05029 - 0.0001537 * t);
  }
  
  
  /**
   * Calculate the apparent longitude of the sun
   *
   * @param t Number of Julian centuries since J2000.0
   * @return The apparent longitude of the Sun in degrees
   */
  private static double sunsApparentLongitude(double t) {
    return sunsTrueLongitude(t) - 0.00569 - 0.00478
        * sin(rad(125.04 - 1934.136 * t));
  }
  
  
  /**
   * Calculate the true longitude of the sun
   *
   * @param t Number of Julian centuries since J2000.0
   * @return The Sun's true longitude in degrees
   */
  private static double sunsTrueLongitude(double t) {
    return geomMeanLongSun(t) + equationOfCentreForSun(t);
  }
  
  
  /**
   * Calculate the equation of centre for the Sun
   *
   * @param t Number of Julian centuries since J2000.0
   * @return The equation of centre for the Sun in degrees
   */
  private static double equationOfCentreForSun(double t) {
    double m = geometricMeanAnomalyOfSun(t);

    return sin(rad(m)) * (1.914602 - t * (0.004817 + 0.000014 * t))
        + sin(2 * rad(m)) * (0.019993 - 0.000101 * t) + sin(3 * rad(m))
        * 0.000289;
  }
  
  
  /**
   * Calculate rad(x) 
   * 
   * @param x
   * @return
   * @since 0.1
   */
  private static double rad(double x) {
    return Math.toRadians(x);
  }
  
  
  /**
   * Calculate deg(x) 
   * 
   * @param x
   * @return
   * @since 0.1
   */
  private static double deg(double x) {
    return Math.toDegrees(x);
  }
  
  
  /**
   * Calculate sin(x) 
   * 
   * @param x
   * @return
   * @since 0.1
   */
  private static double sin(double x) {
    return Math.sin(x);
  }
  
  
  /**
   * Calculate cos(x) 
   * 
   * @param x
   * @return
   * @since 0.1
   */
  private static double cos(double x) {
    return Math.cos(x);
  }
  
  
  /**
   * Calculate tan(x) 
   * 
   * @param x
   * @return
   * @since 0.1
   */
  private static double tan(double x) {
    return Math.tan(x);
  }
  
  
  /**
   * Calculate asin(x) 
   * 
   * @param x
   * @return
   * @since 0.1
   */
  private static double asin(double x) {
    return Math.asin(x);
  }
  
  
  /**
   * Calculate acos(x) 
   * 
   * @param x
   * @return
   * @since 0.1
   */
  private static double acos(double x) {
    return Math.acos(x);
  }
}
